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Abstract 


The viability of a new catchment-based land surface model (LSM) 
developed for use with general circulation models is demonstrated. First, 
simple empirical functions — tractable enough for operational use in the 
LSM — are established that faithfully capture the control of topogra- 
phy on the subgrid variability of soil moisture and the surface water 
budget, as predicted by theory. Next, the full LSM is evaluated offline. 
Using forcing and validation datasets developed for PILPS Phase 2c, 
the minimally calibrated model is shown to reproduce observed evapo- 
ration and runoff fluxes successfully in the Red- Arkansas River Basin. A 
complementary idealized study that employs the range of topographic 
variability seen over North America demonstrates that the simulated 
surface water budget does vary strongly with topography, which can, 
by itself, induce variations in annual evaporation as high as 20%. 
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1 Introduction 


The motivation behind the land surface model (LSM) presented in the com- 
panion paper by Koster et al. (this issue, referred to as Part 1) is the general 
lack of development of runoff formulations compared to sophisticated state-of- 
the-art evaporation formulations [Koster and Milly , 1997]. In particular, most 
LSMs for general circulation models (GCMs) strongly emphasize point pro- 
cesses (e.g. canopy structure and resistances to evaporation) and assume soil 
moisture to be spatially uniform at the GCM scale, the soil being described 
as a stack of vertical layers. The transfer of moisture between these layers, 
which has a crucial impact on both baseflow and overland flow production, is 
then computed one-dimensionally along the vertical. This point-process focus 
neglects the critical effects of the horizontal variability of soil moisture, which 
is largely imposed by precipitation heterogeneity and the lateral redistribution 
of moisture following topography. 

Part 1 proposes a catchment-based LSM for GCMs that accounts explic- 
itly for the control of topography on the small-scale variability of soil moisture 
and its effects on evaporation and runoff. This LSM, which follows on re- 
cent advances in macroscale catchment modeling [Famiglietti and Wood, 1994; 
Stieglitz et al., 1997], recognizes the hydrological catchment as the fundamen- 
tal land-surface unit. In each catchment, the LSM combines a catchment 
model describing the redistribution of water according to topography under 
the effect of gravity with more traditional parameterizations for energy budget 
processes, namely, vegetation and evaporation parameterizations from the Mo- 
saic LSM [Koster and Suarez, 1992; Koster and Suarez, 1996], a snow model 
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from Lynch-Stieglitz [1994], and a linear soil-heat diffusion scheme. 

In the present paper, we examine two critical issues related to the LSM’s 
operational use. First, in Section 2, we determine the extent to which the com- 
plex relationships between topography and the surface water budget described 
in Part 1 can be characterized in terms of empirical functions that are compu- 
tationally manageable. These functions describe, in particular, the impact of 
topography on three processes: (1) baseflow, (2) the water transfers between 
the root zone moisture and the water table, and (3) the lateral distribution 
of the root zone wetness, crucial to the spatial partitioning of the catchment 
into three fractions, each representing a distinct hydrological regime and each 
thereby employing a different treatment of runoff and evaporation. In evaluat- 
ing these basin-specific empirical functions, we also examine (Section 2.6) the 
suitability of automated parameter estimation algorithms, which are critical 
to a global application of the LSM. 

The second issue examined in this paper is the success of the LSM in 
reproducing observed surface fluxes given realistic atmospheric forcing. Such 
an evaluation is, of course, an essential pre-requisite to the full coupling of 
the LSM to an atmospheric GCM. Section 3 presents the framework for an 
offline evaluation study in the Red-Arkansas river basin, adapted from the 
framework used in the PILPS 2C model intercomparison study [ Wood et a L, 
1998]. The fluxes simulated in this environment with the catchment model 
are then compared with observational data in section 4. This section also 
highlights some unique features of the LSM’s behavior present in these results, 
and it presents an additional idealized calculation that very clearly isolates the 
impact of topographical variations on the surface fluxes. 
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Together, Sections 2 and 4 serve to demonstrate that the catchment 
LSM described in Part 1 is more than a set of unconventional but untested 
ideas — the version of the catchment LSM used herein efficiently captures the 
effects of topography on the surface water budget and realistically transforms 
atmospheric forcing variables into surface fluxes. This is emphasized again in 
section 5, which provides some additional discussion. 

2 Topography-related Par ameterizat ions 

In this section, the strategy employed to parameterize the three topography- 
related processes as a function of the catchment moisture and topography is 
first detailed in the case of a specific catchment, the Sleepers River catchment 
of Vermont. This strategy is then evaluated across about 5000 catchments 
spanning North-America. 

2.1 The Sleepers River Catchment 

The hydrological catchment examined first is sub-catchment W-3 (8.4 km 2 ) of 
the Sleepers River Watershed (111 km 2 ), located in the highlands of Vermont 
and monitored since 1957 by the Agricultural Research Service (ARS). More 
details can be found in Stieglitz et a 1. [1997] and Anderson [1977], In partic- 
ular, the topography of this catchment is described by a 30 m x 30 m DEM. 
The topographic index distribution was computed at this resolution [Stieglitz 
et a 1., 1997], and the first three moments are given in Table 1. 

This table also gives the values used to describe soil properties. The 
saturated conductivity at the surface K s (d = 0) was established by Stieglitz 
et al. (1997) for the Sleepers River catchment, based on the values proposed 
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by Beven [1982] for different vegetation cover and soil texture. The param- 
eter v describing the decay of K s with depth was also evaluated by Stieglitz 
et al (1997) for the Sleepers River catchment. The soil matrix potential at 
saturation ip s and the soil parameter b have standard values from Cosby et al. 
[1984], 

2.2 Background: Topography and Water Table Depth 

The TOPMODEL framework gives an analytical relationship between the local 
topographic index x (defined as the ratio of contributing area to local slope; 
see Beven and Kirkby [1979]) and the local water table depth d: 

d = d — ~ (x ~x) (1) 

where x is the mean of the topographic index in the catchment, d is the mean 
water table depth in the catchment, and v describes the exponential decay 
of the saturated hydraulic conductivity K s with depth. Following Sivapalan 
et al. [1987], the actual discrete distribution of the topographic index can be 
conveniently represented with an three-parameter gamma distribution: 

h{x) = Ujjx 0^) exp ^’ x ’ <6>0 ( 2 ) 

The three parameters are derived from the mean (x), variance(cr^), and skew- 
ness {y x ) of the actual distribution of the topographic index, so that the gamma 
distribution has the same first three moments as the actual distribution. This 
imposes that: 

X = 7x^x/2 (3) 
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The resulting three-parameter gamma distribution is displayed in Figure la 
for the Sleepers River catchment. Being analytical, this idealized distribution 
has the advantage of being much easier to manipulate than the actual discrete 
distribution. In particular, the combination of (1) and (2) through a variable 
change leads to an equation for the distribution of the water table depth, given 
its mean d: 


fo{d) = 
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with xd = — xf f an d Vd — d+(x — This distribution, which is reversed 
compared to the gamma distribution of x, is displayed in Figure lb for four 
different values of d. This figure shows that the mean water table depth d 
does not change the shape of the distribution, but just induces a shift along 
the x-axis. Note that the TOPMODEL framework allows “negative” water 
table depths; these correspond to a water table virtually above the surface, i.e. 
to local saturation. 


2.3 Mean Water Depth and Baseflow 

The catchment deficit Mo, an important prognostic variable of the LSM, is 
based on the equilibrium profile w(z) defined in Part 1 (equation 1.4). It 
results from the integration of 1 — w(z), first vertically from the local water 
table depth d to the ground surface, then laterally, over the distribution of 
d across the catchment. The resulting one-to-one relationship between the 
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catchment deficit and the mean water table depth d, which directly controls 
the baseflow production (equation 16, Part 1), can not be written analytically, 
but it can be numerically evaluated. The solid curve in Figure 3a was hence 
constructed using 50 values of d spaced every 10 cm, from d=- 0.9 m to d = 4 m. 

In the catchment-based LSM, this theoretical relationship is approxi- 
mated with a simple analytical function: 

\ 

M d = A (d + Bf (7) 

where the parameters A and B implicitly describe the influence of topography 
on the position of the water table and thus on baseflow. To determine the 
values of A and B in (7), we fit this equation to the theoretical values of Mp 
computed for two values of d, namely d = 0 m and d= 2.5 m. Figure 3a compares 
the resulting approximate function to the theoretical relationship between Mp 
and d. The accuracy of the approximate function is very good for values of 
d < 2.5 m. Note that, in any case, Mp has an absolute upper limit, to reflect 
the existence of bedrock (Part 1); at some point, the growing errors suggested 
at higher water table depths become irrelevant. 

As discussed in Part 1, the inversion of this approximate function pro- 
vides an approximate mean water table depth as a function of the prognostic 
variable Mp. This depth can then be used to approximate the baseflow G 
(using equation 16, Part 1). Figure 3b shows a comparison of this baseflow 
G computed from Mp with that computed directly from d. For catchment 
deficits greater than 100 mm (d > 0.9 m), both the approximate and theoret- 
ical baseflow are smaller than 0.01 mm/h, and the corresponding inaccuracies 
shall have a negligible incidence on the baseflow modeled in the LSM. For 
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smaller catchment deficits, Alp can be successfully used to estimate G. The 
largest differences occur at very low catchment deficit (< 10 mm). Both the 
approximate and theoretical baseflow for this situation are, however, very high 
(> 10 mm/h), so that the catchment deficit is quickly brought back to larger 
values. In any case, the TOPMODEL assumptions leading to the baseflow 
equation are less appropriate at such small values of catchment deficit. 

2.4 Spatial Distribution of the Root Zone Wetness 

2.4.1 Towards an Approximate Distribution 

The spatial distribution of the equilibrium root zone wetness, 8, results from 
the vertical integration across the root zone of the equilibrium soil moisture 
profile, w(z), at each point in the catchment. As discussed in Part 1, this 
distribution of 8 is eventually shifted based on the root zone excess variable 
M r z to account for departures from equilibrium. The resulting distribution 
is then used to separate the unit catchment into three subunits representing 
different hydrological regimes. 

At any point in the catchment, the value of 8 depends non-linearly on 
the position of the local water table depth d with respect to both the ground 
surface and the chosen root zone depth d rz . Under this framework, a tractable 
analytical form for the probability distribution of 8 across the catchment can- 
not be derived from (6). However, this distribution is uniquely related to d 
(or Alp) for a given catchment topography and d TZ , and it can be constructed 
numerically, as shown in Figure 2a for four different values of d in the Sleepers 
River catchment. 

Our strategy is to define an analytical approximation to the distribution 
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of 9 , with parameters dependent on the catchment deficit on one hand, and 
the first three moments of the topographic index on the other. The selected 
distribution is based on a two-parameter exponential distribution: 

f(9) = a 2 (9 - 9 0 ) exp -0 ^ - ^, a, 9 0 > 0 (8) 

where 9o is the minimum value of the distribution, and where a is a shape 

\ 

parameter. The comparison of Figure 2a and 2b indicates that the theoret- 
ical root zone wetness distribution is well approximated by a two-parameter 
exponential distribution, with decreasing a as the catchment deficit decreases. 

At very low catchment deficit, however, the theoretical root zone wetness 
distribution is strongly skewed to the left, and cannot be properly approxi- 
mated by the two-parameter exponential distribution, which is always skewed 
to the right. This deficiency is an important concern for the representation of 
the saturated fraction, A sat , which would be poorly represented by the inte- 
gration of the two-parameter exponential distribution for 9 > 1. Because of 
the crucial role of variable A sat in the present LSM, it was decided to define 
it as a third parameter of the root zone wetness distribution. The complete 
form for the approximate distribution of the root zone wetness at equilibrium 
can thus be written: 

/prox(^) = 0 , 0 < 9 < 9q 

< /prox(0) = A x (9 - 9 0 ) exp~ Q ^ -fl °), a,9 0 > 0; 9 0 < 9 < 1 (9) 

, /prox(l) = A sat 

The term A x is a scaling factor allowing fj f pTOX {9)d9 = 1, and is a trivial 
function of A, at , a and 9 0 . In particular when A, at = 0, A x = 1 and the 
distribution / prox is identical to the the two-parameter exponential distribution. 


10 


2.4.2 Characterization of the Parameters 


The parameterization of the root zone wetness distribution as a function of 
topography and catchment deficit is therefore reduced to the parameterization 
of the three parameters 0q, A sat and a as a function of the topography and 
catchment deficit. The parameters 9o and T sat are easily extracted from the- 
oretical distributions of the root zone wetness, as computed numerically over 
a broad range of water table depths. Figure 4a and 4b show the dependence 
of 0 O and on the catchment deficit M D> in the Sleepers River catchment. 
The approximate curves on these figures are described in both cases by the 
following function: 


Y(M d ) = + (1 - Y x ) — — i±-dy^. 


( 10 ) 


where Y designates either 0 O or A sat . The parameter Tx, describes the position 
of the horizontal asymptote when M D -» oc. In a completely dry catchment, 
there is no saturated fraction, so that 1^=0 when Y = A sat . Because there 
always remains some residual moisture in the soil, Too is, in contrast, positive 
when T = 0 q; in this case, we take T » to be the minimum root zone wetness in 
the distribution constructed for 3=10 m. For both A sat and 0 O , the other three 
parameters A, B, and C are established by Gaussian elimination. Figure 4a 
and 4b show that the resulting approximate functions describe extremely well 
the variations of 0 O and A sat with M D . 

As mentioned in section 2.4.1, the two-parameter exponential distribu- 
tion is an excellent approximation for the actual distribution of the root zone 
wetness 0 for larger values of the catchment deficit. The parameter a can then 
be estimated from 0 mo d) the location of the first mode (densest value) of the 
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theoretical distribution: 


a = 


1 

<?0 0mod 


( 11 ) 


On the other hand, it is not possible to derive a meaningful value for a at low 
catchment deficits, when the three-parameter exponential distribution gives a 
poor representation of the theoretical distribution of 6. Because this distri- 
bution gets flatter, however, as the catchment deficit decreases, we choose to 
define a as an increasing linear function of the catchment deficit (Figure 4c), 
defined by: 


ct — A a Md + B a 


( 12 ) 


The basin-specific parameters A a and B a are empirically fitted, as detailed in 
Appendix A. 3 

Figure 5 compares the resulting distribution / prox (9) to the theoretical 
distribution of the root zone wetness in the Sleepers River catchment. The 
agreement is excellent at high catchment deficit (d = 2 and 3 m), and it 
remains quite adequate even at very low catchment deficit. 


2.5 Transfers between Soil Moisture Prognostic Vari- 
ables 

As described in Part 1, the transfer of moisture between the root zone excess, 
A/rz, and the catchment deficit, Md, requires a timescale of moisture transfer, 
Ti. An empirical equation for T\ as a function of Mrz and Md is obtained 
through offline distributed calculations of subsurface catchment moisture dy- 
namics. In each offline calculation, a highly discretized vertical soil column is 
initialized with a water table depth, an associated equilibrium soil moisture 
profile, and a root zone excess. The soil moisture profile then evolves to a 
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new state over the course of 24 hours using a discretization of the Richards 
equation. For the given one-dimensional initialization, the decrease in the root 
zone excess over the course of a day is thus precisely determined. 

Using the same initial root zone excess, the one-dimensional calculation 
is repeated for a number of different water table depths. Each calculation 
produces a daily loss associated with a specific depth. Now consider that 
a given catchment deficit, iV/ D , and the topographical characteristics of the 
catchment examined imply a unique distribution of water depth within the 
catchment. The daily losses generated in the one-dimensional calculations 
are integrated over this distribution, resulting in a mean loss per day for the 
catchment as a whole, a function of M D and M RZ . n is taken to be the 
timescale consistent with this daily loss. 

The empirical equation that relates T\ to Mr z and M D is: 

n = e (an+fcriiV/ D ), (13) 

where a r \ and b T i depend on the ratio of Mrz to the total water holding 
capacity of the root zone. If this ratio exceeds .01, then a Tl and 6 ri are 
assigned one set of values (determined by fitting to the results of the distributed 
calculations), and if the ratio is less than -.01, then a Tl and 6 Tl are assigned 
a different set of values. Linear interpolation between these two sets provides 
a T i and b rX for intermediate values of the ratio. Recall that tj applies to a 
downward flux (out of the root zone) when Mrz is positive and to a negative 
flux (into the root zone) when Mrz is negative. 

The effectiveness of (13) is illustrated in Figure 6. The solid lines rep- 
resent the timescales derived directly from the distributed calculations (for 
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the case of the Sleepers River catchment), and the dotted lines represent the 
timescales derived with (13). The agreement is generally strong, particularly 
for lower values of iV/ D . As expected, the timescale increases with increasing 
catchment deficit and decreasing root zone excess. 

As discussed in Part 1, the calculation of r 2 , the timescale of moisture 
transfer between the surface excess and the root zone excess, is determined via 
a one-dimensional calculation. The fitted empirical function used to compute 
r 2 has the simple form 

Q t 2 

7 2 = (M rz + 6 T2 M se ) 3 ’ (14) 

where a r2 and 6 r2 are fitted parameters. 

2.6 Adequacy of the Method in North America 

The above section describes the methods used to estimate the model parame- 
ters as a function of the catchment deficit, the influence of topography being 
taken into account implicitly by the parameter themselves. Given the number 
of catchments considered in a global application of this LSM, automation of 
these estimation methods is, of course, essential. In the present section, we 
address the accuracy of our empirical functions over a wide range of catch- 
ments, subject to the application of automated parameter estimation. Note 
that in this section, accuracy refers simply to agreement with the results of 
the more complex, physically-based calculations and not to agreement with 
observations. 

In each of the 5020 catchments covering North- America (Part 1), the 
topographic index has been computed for every 30-arc-second (~ 1 km) pixel, 
based on the topographic information provided at this resolution by the USGS 
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EROS Data Center. Each resulting distribution of the topographic index 
was analyzed in terms of its first three moments. Some combinations of 
moments (e.g. those with negative skew) were not directly amenable to a 
three-parameter gamma distribution and required a preliminary treatment, as 
described in appendix A.l. We also used a scale correction on the catchment 
means of the topographic index (appendix A. 2), based on the recent work of 
Wolock and McCabe [1999]. In all catchments, the required parameters were 
finally estimated by fitting (7), (10) and (12) to the theoretical quantities. The 
automation of this procedure is detailed in appendix A. 3. 

At this stage, an important concern for use in the LSM is that the re- 
sulting approximate baseflow and root zone wetness distribution are accurate 
enough, at least relative to the theoretical quantities. 

2.6.1 Accuracy of Approximate Baseflow 

The accuracy of the approximated baseflow was quantified by computing the 
Root Mean Square Error (RMSE) between it and its theoretical value (i.e. 
computed directly from d). In every catchment, the largest contribution to 
this RMSE is from the inaccuracies at very small catchment deficit. Because 
such small Mp are rare in the LSM, the RMSE was only computed for values of 
baseflow corresponding to M D > 2 mm. This RMSE was computed in all 5020 
catchments in North-America, and was found, in 83% of them, to be smaller 
than 0.15 mm/h, the value computed in the Sleepers River catchment, where 
baseflow is very successfully approximated (Figure 3b). The largest RMSE 
found across North-America is 1.13 mm/h, and Figure 7 shows a comparison 
of the approximate and theoretical baseflow in the corresponding catchment. 
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Despite large errors at very small catchment deficits, the overall behavior of 
the approximate baseflow in this “worst case” still closely follows the theoret- 
ical behavior, which is ruled by topography. This study shows the baseflow 
RMSE to be reasonably small in the 5020 catchments, which demonstrates the 
adequacy of the automated parameter estimation procedure. 

2.6.2 Accuracy of Approximate Root Zone Distribution 

We focused here on the accuracy of two quantities derived from the root zone 
wetness distribution and used in the LSM. The first one is the saturated frac- 
tion A sat itself, which is used for the areal partitioning of each catchment. 
The second quantity is the mean root zone wetness in the unsaturated but 
“unstressed transpiration” fraction, defined as: 

9tT= L ^/pro *(6)d0 / f fprox(0)d0 (15) 

where 0 w nt is the wetness at the wilting point. Figure 4b shows, for the Sleepers 
River catchment, an excellent agreement between the theoretical and approx- 
imate values of the saturated fraction A sa t. As for baseflow, the largest errors 
in A sat occur at very low catchment deficit. Again, these low values being very 
rare in the LSM framework, the RMSEs of both A sat and 9^ were computed 
in all 5020 catchments for values of M D > 2 mm. 

This analysis shows a satisfying accuracy of .4 sat in all catchments: 50% 
of the catchments have a RMSE smaller than the RMSE of A sat in the Sleepers 
River catchment (RMS=0.032) where A sat was particularly well approximated 
(Figure 4b). The largest RMSE for A sat is 0.093, and Figure 8a compares the 
approximate and theoretical saturated fraction in the corresponding “worst 
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case” catchment. The agreement between the two curves is excellent, despite 
errors at low values of catchment deficit. In particular, the RMSE drops to 
0.010 when the values corresponding to M D < 15 mm are excluded from the 
RMSE computation. 

Figure 8b shows a very satisfying agreement between the approximate 
and theoretical values of 0 tr in a catchment with RMSE = 0.018. This RMSE 
is exceeded in only 17 % of the catchments, which demonstrates the general 
accuracy of the approximate (fc. Finally, Figure 8c compares the approximate 
and theoretical 0 tr in the catchment with the highest RMSE, where the un- 
derestimation of the approximate $i T is clearly significant. Still, this degree of 
inconsistency is the exception rather than the rule, and it is not, in any case, 
large enough to have a first order effect on the results. 

2.6.3 Accuracy of Approximate Transfers 

At this time, the automated procedure used to extract the parameter values 
used in the calculation of r u via (13), has been tested on a much more limited 
scale. For a variety of randomly sampled topographic conditions, the procedure 
generates parameters that reproduce the timescales obtained through highly 
distributed calculations. The success in Figure 6 is typical of that obtained in 
these tests. 

As for the calculation of r 2 , the present version of the model uses the 
same equation and parameters in every catchment; we assume that moisture 
transfer between the surface layer and the root zone is not, to first order, 
strongly affected by the topography. Thus, a test of the automated procedure 
for estimating the parameters in (14) is not performed. 
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3 The Offline Testing Framework 


The accuracy of an LSM can only be determined through a detailed com- 
parison of its products with observations. However, the testing of an LSM 
coupled to an atmospheric model can be very difficult if the forcing by the at- 
mospheric model is in error. For instance, GCM-simulated precipitation rates 
are notoriously inaccurate, and these errors would probably dominate the be- 
havior of the LSM. This explains the widespread use in validation studies of 
“offline” simulations, for which the atmospheric forcing is derived from ob- 
servations [Henderson-Sellers et al., 1996; Shao et a l, 1994; Dirmeyer et al., 
1999], 

3.1 Description of the Offline Framework 

The offline framework developed over the Red-Arkansas River Basin for the 
Project for Intercomparison of Land-surface Parameterization Schemes (PILPS) 
Phase 2c [Wood et al., 1998] was adapted for the catchment-based LSM. Fig- 
ure 9 compares the boundaries of the Red-Arkansas River Basin in PILPS2c 
and in our study. The latter are much more realistic, since our unit catchments 
are extracted from a DEM. In particular, the total area is much closer to the 
actual area proposed by the USGS (566,196 km 2 ) in the LSM (area = 571,890 
km 2 ) than in the original PILPS2c framework (area = 608,211 km 2 ). The fig- 
ure also illustrates the gain in average resolution, as the average area of the unit 
catchments (4540 km 2 ) is about half the average size of the l°x l°PILPS2c 
grid-cells (9970 km 2 ). 

Most forcing variables were provided at the l°x 1° resolution over the 
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PILPS2c domain (Figure 9). The meteorological forcing, available at a 1-hour 
resolution for ten years (1979-1988), consists of air temperature and humidity 
at two meters, surface wind and pressure, precipitation, and downwelling solar 
and longwave radiation. These variables were interpolated from observational 
point data from the National Climatic Data Center (NCDC). The vegetation 
properties (leaf area index, greenness fraction, roughness length and snow-free 
albedo) are derived from the ISLSCP data set for 1987-1988 at the l°x 1° res- 
olution [Sellers et a/., 1996]. 

The averages of the 1987 and 1988 monthly values are used; thus, sea- 
sonal variations in vegetation properties are included but interannual variations 
are neglected. Porosity and depth of the soil profile have been interpolated 
from a 1-km resolution soil data base [STATSGO, 1994] to the l°x 1° resolu- 
tion. Because of the mismatch between the catchment space and the regular 
lxl grid, the above forcing fields have been interpolated to the catchment 
space. To this end, the intersections between the catchments (defined as poly- 
gons) and the grid cells have been computed by triangulation, allowing the 
simple areal weighting interpolation explained in Part L 

The vegetation type in each catchment w'as assumed to be one of the 
eight surface types used in the Mosaic LSM [Koster and Suarez, 1996]. The 
type represented the most within the catchment’s area, as determined from 
vegetation maps on the ISLSCP dataset, was assumed to cover the catchment 
uniformly. As for the soil properties not mentioned above, they were taken 
from Table 1 in all catchments. 

The initial conditions for the LSM were obtained by K spinning up” the 
model for 1 year (with the forcing from 1979) before running it for the ten 
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years. There is therefore no strong equilibrium constraint on the LSM, and 
this enhances the comparability between modeled and observational quanti- 
ties [Wood et ah, 1998]. 

3.2 Validation Data 

Relative to most other LSM validation frameworks, the “Red-Arkansas” frame- 
work is spatially extensive and features a lengthy validation period. This 
multi-year dataset allows us to assess the interannual variability of the fluxes 
generated by the LSM. Another interesting feature is the quality of the vali- 
dation data, which include (1) daily naturalized streamflows [ Lohmann et a l, 
1998] at the two downstream gauging stations (Little Rock and Shreveport, see 
Figure 9), and (2) monthly total evaporation over the entire domain, estimated 
from 1980 to 1986 through an atmospheric budget analysis using radiosonde 
data [Liang et a 1., 1998]. Over the 80-86 period, the accumulated sum of these 
two validation data sets balances very closely the accumulated prescribed pre- 
cipitation, which demonstrates, under the assumption of minimal long-term 
soil moisture change, the consistency between these two independent data sets 
and the forcing data. 

A routing scheme was used in PILPS2c [Lohmann et al., 1998] to simu- 
late streamflow from modeled runoff. This scheme, designed to route regularly 
gridded runoff, cannot be easily applied to the catchment-based LSM’s prod- 
ucts. Because the delay of the routed relative to the non-routed runoff peaks 
is only 3 to 5 days (Lohmann, personal communication), we can avoid the 
routing issue by focusing on the monthly to annual timescales. 
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4 Numerical Results 

4.1 Simulation of Fluxes in the Red- Arkansas 

Our comparisons with observations in the Red- Arkansas basin focus mainly on 
the seasonal cycles of areally-averaged evaporation and runoff. Figures lOab 
show that the model reproduces well the mean seasonal cycles of these two 
quantities. Runoff is overestimated slightly in the fall, and evaporation is 
underestimated slightly in the summer, but these errors are probably of the 
same order as the observational error itself. We must emphasize, however, that 
the Red-Arkansas results were used to a certain extent in the development of 
this model; thus, the agreement with observed fluxes does not constitute a 
model validation. We present Figures lOab mostly to show that the model, 
with its current design, can reproduce the observed fluxes without an arbitrary 
tuning of nebulous model parameters. Indeed, the only parameter that was 
“tuned” in the exercise was the soil depth associated with the surface excess 
variable. The chosen value for this depth, 5 cm, is consistent with surface layer 
thicknesses employed by many modeling groups. 

Figures lOcd show that the model reproduces the interannual variations 
of the mean runoff and evaporation cycles as well. Normalized fluxes are 
shown. The monthly dynamics of total runoff, which are mainly driven by the 
precipitation signal, are realistically simulated: high and low values occur in 
the correct months, and the amplitude of the seasonal cycle is realistic, de- 
spite a tendency to overestimate runoff in dry periods (e.g. the fall of 1980 
and the summer of 1981) and conversely to underestimate it during the high 
runoff periods (e.g. winter, 1984/85). Some of the differences between the 
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simulated runoff and its observational surrogate (streamflow divided by con- 
tributing area) may be related to the lack of runoff routing in the model. The 
simulated runoff is necessarily more variable than the routed runoff, being well 
correlated to precipitation. (For example, both precipitation and simulated 
runoff drop sharply in January, 1985.) The high peak for the observed nor- 
malized runoff in January 1986 reflects the very small precipitation in that 
month; the absolute runoff error in that month is thus quite small. 

4.2 Unique Aspects of the Model Behavior 

Aside from validation, the Red-Arkansas framework allows us to demonstrate 
original features of the catchment-based LSM’s behavior. Figure 11 illustrates 
the influence of the meteorological conditions on the catchment partitioning, 
which underlies the new LSM. The 1980-86 mean precipitation in catchment 
A, which is in the eastern part of the river basin, is 1246 mm/y, whereas 
that in catchment B, in the western part, is 635 mm/y. The figure shows 
that this difference in precipitation leads to different seasonal dynamics of 
the wilting, transpiration, and saturated areal fractions. In both catchments, 
the saturated fraction is higher in winter (when precipitation exceeds evap- 
oration) and decreases in summer (when evaporation exceeds precipitation), 
with a corresponding summer increase in the wilting fraction. In the wet 
catchment, however, the saturated fraction can reach 13% in March, and the 
wilting fraction appears only during summertime. In contrast, the wilting frac- 
tion is greater than 30% all year long in the drier catchment, and the saturated 
fraction there is always very small. 
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4.3 An Additional, Idealized Calculation 


A clearer demonstration of topographic effects is afforded by the topograph- 
ical description of North American catchments presented in Section 2. In 
an idealized experiment, we force all 5020 catchments with the same single 
timeseries of meteorological forcing. All catchments are also given the same 
surface parameters, such as vegetation type (grassland), leaf area index, and 
roughness length. Thus, the only differences between the catchments are their 
topographical statistics and their associated fitted functions for hydrological 
regime separation, baseflow generation, and subsurface moisture transport. By 
comparing the variations in the fluxes generated by the different catchments, 
we can get a sense for how important the modeled topographic effects are, 
given the range of topographical characteristics likely to be encountered in a 
global simulation. 

The simulation is broken up into three time segments. In the first seg- 
ment, a year of meteorological forcing is applied repetitively until all catchment 
moisture variables come to equilibrium, i.e., until the variables attain the same 
values at the end of the year as they have at the beginning. The meteorological 
forcing data used for this year is the 1987 forcing data for a specific catchment 
in the Red-Arkansas basin, as derived from the ISLSCP CD-ROM1 [Sellers c t 
al., 1996]. The same catchment is used to define the surface characteristics of 
all the catchments in the test. The second time segment is a one-year contin- 
uation of the previous simulation, but using the 1988 forcing for the specified 
Red-Arkansas catchment. The second time segment thus represents a dry year 
following a relatively wet year (the annual precipitation was 739 mm and 353 
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mm for 1987 and 1988, respectively). In the third segment, the 1988 forcing 
data is applied repetitively to bring the catchment moisture variables to a new 
equilibrium. 

Annual values of evaporation and runoff from the three time segments 
are normalized by the corresponding annual precipitation and plotted against 
certain topographical statistics in Figure 12. Again, the catchments differ 
only in their topographical description; thus, if a similar experiment were per- 
formed with a standard one-dimensional SVAT model, all catchments would 
produce precisely the same runoff and evaporation values. Figure 12 shows 
that differences in topography can lead to significant differences in runoff (up 
to about 20%). Differences in evaporation can be as high as 5% for the equili- 
brated years and 20% for the transition year (the second time segment). Note 
that evaporation exceeds precipitation during this second year; moisture stored 
during the relatively wet 1987 is evaporated during 1988. 

The topographic statistics chosen for plotting are those that correlate 
best with runoff and evaporation. (The presence of some scatter, of course, 
shows that other statistics can also be important.) Runoff correlates best 
with the skewness of the topographic index; for this imposed climate, overland 
flow dominates over baseflow, and the skewness helps determine the saturated 
fraction of the catchment, which contributes to overland flow. Evaporation, on 
the other hand, correlates best with the standard deviation of the topographic 
index, presumably because the standard deviation contributes most to the 
shape of the root zone moisture distribution at the onset of wilting, which in 
turn determines the dynamics of the wilting fraction, A w n t . 

Figure 12 clearly demonstrates the impact of topographic effects, as mod- 
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eled, on the annual mean surface water balance. The next logical step, of 
course, is to validate these topographic effects against observations. The iso- 
lation of the topographic effects from all other effects (meteorology, surface 
characteristics, etc.) in the observational record will be addressed in a future 
study. 

5 Conclusions and Discussions 

Many studies have demonstrated the impact of small scale spatial variability in 
hydrological processes on the average interactions between the land surface and 
the atmosphere [e.g. Avissar and Pielke, 1989, Johnson et a 1, 1993, Ducharne 
et a 1, 1998]. This impact limits the ability of point-process models to represent 
large scale hydrological processes realistically. These considerations motivated 
the development of a new land surface model for GCMs, the basic framework 
for which is provided in Part 1. This new model calls for the disaggregation of 
the land surface into a mosaic of hydrological catchments, determined through 
the analysis of a high resolution DEM. In each catchment, a catchment model 
describes the redistribution of moisture according to topography under the 
effect of gravity, based on TOPMODEL [Beven and Kirkby, 1979; Sivapalan 
et al., 1987]. The resulting distribution of the root zone moisture is then used 
to subdivide the catchment into three fractions, each representing a distinct 
moisture regime; evaporation and runoff processes are modeled with adapted 
formulations in each fraction, in order to produce a more accurate estimate of 
the catchment-mean rates. 

Three hydrological processes (baseflow, water transfers in the unsatu- 
rated zone, spatial distribution of the root zone moisture) depend on both 
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topography and the bulk catchment moisture in a highly complex way. This 
paper explains in detail how these processes are approximated by simple, 
tractable functions of three moisture prognostic variables. Under this frame- 
work, the influence of topography is implicitly represented by the parameters 
of the fitted functions. These functions are pre-processed in every unit catch- 
ment, and the approximate relationships are used at each time step in the 
model, which strongly enhances the performance of the model in terms of 
memory and CPU requirements. An extensive study amongst the 5020 catch- 
ments covering North-America showed that automated procedures produce 
fitted, empirical functions that capture very well the behavior inherent in the 
more complex, physically-based calculations. Nevertheless, several improve- 
ments can still be brought to the parameter estimation procedure. The fitting 
of the approximate functions to the theoretical relationships could be opti- 
mized, using, for instance, RMSEs as criteria. The chosen form of the fitted 
functions could perhaps also be improved. 

A model using these fitted functions and the other components described 
in Part 1 has been coded and tested for consistency. This model has then been 
run offline in the Red- Arkansas River Basin using the framework developed for 
PILPS2c [Wood et a L, 1998]. The overall performance of the catchment-based 
LSM in this exercise is very satisfactory (Figure 10). Though this does not con- 
stitute model validation (indeed, the structure of the model evolved partly as a 
result of evaluations against these data), the offline simulation does show that 
the model can accurately reproduce, with a reasonable selection of parameter 
values and minimal tuning, the observed seasonal and interannual variations 
in runoff and evaporation. This offline study also demonstrates the dynamic 
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fractioning of the unit hydrological catchments into three sub-fractions (Fig- 
ure 11). The model simulates in particular the saturated fraction in the unit 
catchments, and this has important implications for validation, since it is now 
possible to estimate that quantity through remote-sensing [Gineste et al, 1998; 
Franks et al., 1998]. 

An important implication of this study is that evaporation and runoff in 
a given catchment depend strongly on its topography. This is demonstrated 
most clearly in Figure 12, which shows, for a dry year following a relatively 
wet year, topographically-induced variations in annual evaporation spanning 
20%. Evaporation appears to be controlled mostly by the standard deviation 
of the topographic index, whereas surface runoff is controlled mostly by the 
skew. 

Many directions can still be explored to improve the realism of this 
catchment-based LSM. Because the distribution of the topographic index is 
known to be altered at coarse resolutions (e.g., Wolock and Price, 1994; Zhang 
and Montgomery , 1994; Quinn et al., 1995], algorithms to correct the three mo- 
ments of the topographic index from the effect of the coarse 1-km resolution 
in GTOPO30 (following for instance the work by Wolock and McCabe [1999] 
used in Appendix A. 2), or the use of a finer resolved DEM, are an obvious first 
direction. We suspect that the model’s sensitivity to variations in topography 
may be enhanced when finer-scale topographical data are considered. 

Second, most soil parameters are, in the present incarnation of the model, 
identical in all catchments but should vary geographically. The value of the 
TOPMODEL parameter v, which describes the vertical decay of transmissiv- 
ity, is particulary difficult to prescribe, since it is not directly measurable. As 
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a first approximation, it could be described as a function of active soil depth. 
Third, a spatial distribution of the vegetation cover should eventually be intro- 
duced in each unit catchment, in order to better describe the actual vegetation 
properties from the 8 Mosaic vegetation types. Fourth, as mentioned in Part 1, 
the development of an elevation-based disaggregation scheme for atmospheric 
forcing, such as precipitation rate or temperature, using the topographical in- 
formation available for each catchment, should increase the accuracy of the 
catchment model’s outputs. 

The list of potential improvements to the model, of course, extends even 
beyond those just noted. Despite any current shortcomings, though, the model 
is demonstrated here to be a viable tool for estimating large-scale surface 
energy and water budgets, a tool that could, with a little work on the interface, 
be incorporated into a GCM. The fact that the model is not limited by the 
traditional one-dimensional framework of current LSMs should stand it in good 
stead as evolving GCMs continue to require more realistic representations of 
energy and water budgets at the land surface. 
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A Generation of the Parameters in North Amer- 
ica 

A.l Preliminary Treatment of the Topographic Index 
Distribution 

In 110 catchments amid the 5020 covering North-America (2.2 %), the topo- 
graphic index x has a negative skewness and can therefore not be idealized 
with a three-parameter gamma distribution (positively skewed). A gamma 
distribution fr{x) can however be constructed using the absolute value of the 
skew. A distribution /r„, g (x), with the same first three moments as the topo- 
graphic index x, can then be constructed from fr( x ) by symmetry around x: 
/r neg (a;) = f T ( 2x — x). The mean and variance are those of fr{ x ), and the 
skew is inversed, and therefore negative, as required. 

Another issue is that the gamma distribution is not defined if 7 X = 0 (3, 

4). Since this distribution is numerically approximated [Press et al., 1992], very 
small skewnesses induce numerical problems too, and were therefore modified. 
Only two catchments have a skewness in [-0.005, 0.005]. They were changed 
to ±0.005, depending on the sign of the original skewness, and the skew was 
recomputed, using the variance, to define x from (3). 

Finally, the shape of the gamma distribution for <f> < 1 (or skewness 
> 2), with a vertical asymptote at x = y, (Figure 13), is not adequate for 
the numerical construction of neither the root zone wetness distribution nor 
the relationship between d and M D (7). Because this shape seems, in any 
case, highly unrealistic for the topographic index, skewnesses were reduced 
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to values smaller than 2 (in absolute value). To prevent numerical problems 

associated with the shape of the root zone wetness in the case of topographic 

index distributions that are too sharp, the following adjustments were made: 

if 7 X > 1.9 then 7* = 1.9 

if 7 X < -1.6 then j x = -1.6 

These adjustments addressed 265 catchments(5.3%) in the case of positive 
skewnesses, and 10 catchments (0.2 %) in the case of negative skewnesses. 

A. 2 Scale Correction 

A number of studies have shown that the distribution of the topographic in- 
dex is strongly affected by DEM resolution [e.g., Wolock and Price, 1994; 
Zhang and Montgomery, 1994; Quinn et a l, 1995]. More recently, Wolock and 
McCabe [1999], based on 50 catchments in the conterminous United States, 
proposed an empirical relationship between the mean of the topographic in- 
dex computed at the 100-meter and 1000-meter (or 30-arc-second) resolutions 
(mean 10 o and meaniooo respectively): 

meanioo = —1-957 + 0.961 mean 10 oo (16) 

This linear regression explains 93% of the variance of mean 10 o- Such a correc- 
tion of the distribution of the topographic index, even though limited to the 
mean, was applied to all 5020 catchments in North-America. The sensitivity 
of the LSM to this correction is currently under study. 

A. 3 Automated Parameter Estimation 

The different soil characteristics chosen to parameterize basefiow and the root 
zone wetness distribution in each catchment are, for simplicity, taken from 
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Table 1. More detailed soil maps [ e.g ., STATSGO , 1994] could, of course, be 
used to give the parameters a more realistic geographical distribution. 

Equation 7, which describes an approximate relationship between Md 
and d, can easily be fitted in all 5020 unit catchments in North-America, using 
the same approach used in the Sleepers River catchment (section 2.3). 

The parameters for Agat and 0 O (10, 12) are also determined for each of 
the 5020 catchments in North America. The automated procedure, based on 
Gaussian elimination in equation (10), requires three points ( Md,Y ) for both 
Y = A sat and Y = 9q. In both cases and each catchment, these three points 
must meet two important requirements: they must describe the whole range 
of variation of Y, and the resulting second-degree polynomial in denominator 
in (10) must not have its roots in the range of possible values for Mp. 

The automated procedure to estimate a distinguishes two cases according 
to the value of the skewness -y x . If 7 r > 0.25, then the dependence of a on 
Md is described by a straight line, like in the Sleepers River catchment (12). 
The determination of the two parameters A a and B a requires two points (Md, 
a). The first one is (Mp= 0, a=l), and the second one is given by the highest 
considered catchment deficit (for d= 4 m) and the corresponding value of a, 
using (11). If 7 X < 0.25, then the dependence of a on M D is rather described 
by two joining segments: 

if Mq < 100 mm then a = A\ jV/ d + B\ (17) 

if Mq > 100 mm then a = A 2 Md + B 2 (18) 

The two points used to define A a and B a from (12) are kept to define two of 

the unknowns. The remaining two unknowns are given by the intersection of 
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the two segments, for Mp= 100 mm and the corresponding a. 

These two segments allow us to account for the fact that, for low or 
negative skewnesses, the approximate root zone wetness distribution better 
follows the theoretical distribution if the shape parameter increases at a smaller 
rate for small values of the catchment deficit. 

\ 
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Table 1: 


Sleepers River catchment: description of soil and topography. 



Parameter 

Value 

Unit 

Topographic 

Mean x 

7.36 

- 

index 

Variance a 2 

5.51 

- 


Skewness 7 * 

1.23 

- 

Soil properties 

Saturated hydraulic conductivity K s {d — 0 ) 

2.2 10~ 3 

m.s -1 


Transmittivity decay factor v 

3.26 

m _1 


Saturated soil matrix potential V’s 

-0.281 

m 


Clapp and Hornberger (1978) b 

4 

- 


Root zone depth d rz 

1 

m 


Wetness at the wilting point 0 w n t 

0.26 

- 
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d = d — j (x — x) 

I 



d (m) 


Figure 1: Relationships between (a) the theoretical distribution of the topo- 
graphic index x, and (b) the distribution of depth to the water table d, for 4 
values of the mean water depth, in the case of the Sleepers River. 
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Figure 2: Distribution of the root zone wetness at equilibrium, 9: (a) theoret- 
ical distribution for 4 values of the mean water depth in the Sleepers River, 
(b) two-parameter exponential distribution for 9q = 0 and different values of 
a. 
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(a) (b) 



Mean water table depth (m) M 0 (mm) 

Figure 3: Sleepers River catchment: comparison of the theoretical and approx- 
imated^) catchment deficit Mq (mm) as a function of the mean water table 
depth d (m); (b) baseflow G (mm/h) as a function of the catchment deficit 
(mm). 


(a) (b) 



M d (mm) M d (mm) 

Figure 4: Sleepers River catchment: comparison of the theoretical and approx- 
imated values of a) 9q and b) a\, as a function of the catchment deficit. 
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Density 


d_ 0 m d= 1 m d= 2 m d= 3 m 
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9 e e e 

Figure 5: Root zone_vvetness distribution for four different values of the mean 
water table depth ( d= 0, 1, 2, and 3 m) in the Sleepers River catchment: 
comparison of the theoretical and approximated distributions. 
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Figure 6: Contours showing the timescale T\ (hours) as computed from the 
distributed calculations (solid lines) and the empirical approximation, (6). 
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Saturated fraction A** (-) 


RMSE = 1.13 mm/h 



Figure 7: Accuracy of approximated basefiow : comparison of theoretical and 
approximated basefiow (mm/h), as a function of the catchment deficit M D 
(mm), in the unit catchment with the highest RMSE. 



M 0 (mm) 


M 0 (mm) 


M 0 (mm) 


Figure 8: Accuracy of approximated root zone wetness distribution: compari- 
son of theoretical and approximated quantities, as a function of the catchment 
deficit Md (mm) : (a) saturated fraction A sat (-) in the unit catchment with 
the highest RMSE; (b) 0 lr in a unit catchment with RMSE = 0.018, and (c) 
0 tT in the unit catchment with the highest RMSE. 
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Figure 9: Definition of the “Red-Arkansas” spatial domain in the LSM, with 
the Arkansas river basin (light grey, 85 unit catchments) upstream of Little 
Rock (AK) and the Red River basin (dark grey, 41 unit catchments) upstream 
of Shreveport (LA). The thick regular defines the Red-Arkansas domain in 
PILPS2c, at the l°x 1° resolution (61 grid-cells). 
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Figure 10: Comparison of simulated and observational fluxes in the Red- 
Arkansas Basin in 80-86: mean seasonal cycles of (a) evaporation rate and 
(b) total runoff (mm/d); monthly ratios of (c) runoff to precipitation and (d) 
latent heat flux to net radiation (-). The maximum peaks of observational 
ratios (at 7.46 for LE/Rn, and 1.11 for R/P) are truncated in the figure. 


46 






Areal fractions 


Catchment A Catchment B 



J FMAMJ JASOND J FMAMJ J ASOND 

Figure 11: Red-Arkansas Basin: mean seasonal cycle (80-86) of the three 
catchment’s fractions (saturated, unsaturated and unstressed, and stressed 
fraction) in two unit catchments: A: an eastern and humid catchment, and B: 
a western and dry catchment. 
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Figure 12: Left: Variation of annual runoff (normalized by annual precipitation 
with the skew of the topographic index, for three time segments of an idealized 
simulation (see text). Right: Variation of annual evaporation (normalized by 
annual precipitation) as a function of the standard deviation of the topographic 
index. 
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